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CONFERENCE  PAPER 

Situational  awareness  of  Earth-orbiting  particles  is  highly  important  for  future  human  activities  in  space.  For  optical 
observations  of  debris,  multiple  observations  must  be  combined  in  order  to  determine  the  orbit  of  the  observed 
object.  It  is  generally  uncertain,  however,  whether  two  arbitrary  tracks  are  of  the  same  object,  and  solving  this 
problem  can  be  computationally  intensive.  In  this  paper,  we  propose  a  technique  of  correlating  multiple  optical 
observations  by  means  of  probability  distributions  in  Poincare  orbit  element  space. 

1.  INTRODUCTION 

Situational  awareness  of  Earth-orbiting  particles  such  as  active  satellites  and  space  debris  is  highly  important  for 
future  human  activities  in  space.  Presently,  over  300,000  particles  have  been  estimated  to  exist,  and  over  80,000 
observations  are  made  per  day  [5].  Observations  are  made  either  by  radar  or  optical  sensors.  For  optical 
observations,  which  are  usually  made  for  objects  in  medium  Earth  orbit  (MEO)  and  geostationary  orbit  (GEO),  only 
the  angles  and  angular  rates  of  the  track  can  be  determined.  That  is,  the  range  and  range-rate  remain  largely 
unknown.  Therefore,  in  order  to  determine  the  orbit  of  the  observed  object,  multiple  observations  must  be 
combined.  It  is  generally  uncertain,  however,  whether  two  arbitrary  tracks  are  of  the  same  object,  and  solving  this 
problem  can  be  computationally  intensive.  This  is  the  crux  of  the  too  short  arc  (TSA)  problem.  Milani  et  al.  have 
proposed  a  solution  for  heliocentric  orbits  where  each  track  is  expressed  in  a  4-dimensional  quantity  called  the 
attributable  vector ,  and  by  placing  a  few  physical  constraints,  they  restrict  the  range  and  range-rate  to  a  region 
called  the  admissible  region  [4].  Discretized  points  on  the  admissible  region  are  referred  to  as  Virtual  Debris  (VD) 
particles.  Tommei  et  al.  expanded  this  method  to  Earth  orbiting  objects  [6].  Maruskin  et  al.  introduced  another 
method  that  uses  maps  of  the  admissible  region  in  Delaunay  orbit  element  space  [3]. 

In  this  paper,  we  propose  a  technique  of  correlating  multiple  optical  observations  by  means  of  probability 
distributions  in  Poincare  orbit  element  space.  We  first  define  the  admissible  region  mathematically,  as  well  as 
introduce  other  necessary  concepts  (Section  2).  We  then  explain  our  method  and  how  we  incorporate  observation 
data  (Section  3).  An  admissible  region  for  an  observation  is  mapped  to  the  6-D  Poincare  space,  which  is  discretized 
into  many  hypercubes  or  bins.  At  each  bin,  the  density  of  VD’s  is  determined,  and  its  distribution  over  the  Poincare 
space  can  be  regarded  as  a  probability  density  function  (pdf)  as  to  where  the  observed  object  may  exist.  We 
combine  pdfs  from  multiple  observations  using  Bayes'  theorem.  A  significant  computational  bottleneck  in  this 
proposed  method  arises  from  the  very  large  number  of  VD's  that  must  be  mapped  non-linearly  in  order  to 
completely  represent  the  admissible  region  in  6-D  space.  We  avoid  this  problem  by  approximating  the  admissible 
region  as  a  conglomerate  of  smaller  subsets,  and  linearly  mapping  these  regions.  Finally,  we  discuss  a  MATLAB 
implementation  of  our  method  (Section  4).  We  simulated  correlating  996  error- free  optical  observations  for  8 
objects  in  MEO  and  GEO  over  the  course  of  approximately  24  hours.  All  observations  were  correctly  correlated 
with  no  false  positives.  Even  in  the  presence  of  observation  error,  Monte  Carlo-like  test  results  suggest  that  the 
method  will  perform  well.  We  also  tested  the  code  in  more  difficult  observation  geometries  to  ascertain  its  limits 
(Section  5).  Namely,  we  considered  the  case  where  2  objects  lie  within  a  discretization  unit,  where  1  object  is 
observed  simultaneously  at  2  observatories,  and  where  a  satellite  constellation  is  observed  at  1  observatory. 

2.  BACKGROUND 

In  this  section,  we  introduce  the  mathematical  definition  of  the  attributable  vector  and  the  admissible  region,  and 
both  the  exact  (non-linear)  and  linearized  transformations  from  topocentric  spherical  coordinates  to  Poincare  orbit 
elements.  We  also  discuss  the  topology  of  admissible  regions,  which  are  2-dimensional  manifolds  embedded  in  6- 
dimensional  state  space.  Unless  otherwise  stated,  all  length  units  in  this  paper  are  in  Earth  radii  rE,  time  units  in 
hours,  mass  units  in  kilograms,  and  angles  in  radians. 


Form  Approved 
OMB  No.  0704-0188 


Report  Documentation  Page 


Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 


1.  REPORT  DATE 

SEP  2010 

2.  REPORT  TYPE 

3.  DATES  COVERED 

00-00-2010  to  00-00-2010 

4.  TITLE  AND  SUBTITLE 

Correlation  and  Initial  Orbit  Determination  for  Short- Arc  Optical 
Observations 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

University  of  Colorado  at  Boulder, Dept,  of  Aerospace  Engineering 
Sciences, Boulder, CO, 80309 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS (ES) 

10.  SPONSOR/MONITOR’S  ACRONYM(S) 

11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

2010  Advanced  Maui  Optical  and  Space  Surveillance  Technologies  Conference,  14-17  Sep,  Maui,  HI. 
Sponsored  in  part  by  AFRL.  U.S.  Government  or  Federal  Rights  License 


14.  ABSTRACT 


15.  SUBJECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 
ABSTRACT 

18.  NUMBER 
OF  PAGES 

19a.  NAME  OF 
RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Same  as 
Report  (SAR) 

10 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


2.1  The  Attributable  Vector 

For  optical-only  observations,  which  are  usually  made  for  objects  in  medium  Earth  orbit  (MEO)  and  geostationary 
orbit  (GEO),  only  the  angles  and  angular  rates  of  the  track  can  be  determined  [3].  That  is,  the  range  and  range  rate 
remains  largely  unconstrained,  except  for  a  few  physical  restrictions  which  can  be  used  to  constrain  their  values. 
Thus,  each  track  can  be  mathematically  expressed  in  terms  of  an  attributable  vector  ft  at  epoch  t  of  the  observation 
[6]: 


ft  =  (a,  6 ,  a,  6)  e  [-tt,  n)  x  (- tt/2 ,  tt/2)  x  R2,  ( 1 ) 

where  a  and  S  specify  the  topocentric  angular  position  of  the  debris  particle.  A  discussion  of  how  one  may  estimate 
an  attributable  vector  from  a  given  track  of  data  can  be  found  in  Maruskin,  et  al  [3].  We  use  J2000  as  our  coordinate 
system  so  that  a  is  the  right  ascension  and  S  is  the  declination.  In  addition,  information  regarding  time  and  the 
location  of  the  observer  should  be  stored  for  a  more  complete  description  of  the  track,  leading  to  an  extended  set  X: 

X  =  (ftj0  ,M,®)  (2) 


where  t0  is  the  time  of  the  observation,  h  is  the  altitude  of  location  of  observation,  and  0  and  cp  are  the  angular 
position  of  observation  for  a  geocentric  spherical  coordinate  system.  We  chose  a  coordinate  system  such  that  0  is 
the  latitude  and  (p  is  the  longitude  of  the  observation  point.  In  the  following  discussion,  we  will  ignore  h. 

2.2  The  Admissible  Region 

For  some  attributable  vector  X ,  we  can  take  different  values  of  range  and  range-rate  (p,p)  to  complete  the  topocentric 
coordinates  of  the  particle  and  thus  obtain  different  physical  orbits.  Visually,  we  can  imagine  taking  different  points 
in  the  topocentric  range  /  range-rate  plane.  However,  not  all  of  these  orbits  are  relevant  for  any  given  application. 
Rather,  a  closed  region  of  the  (p,p)  plane  can  be  defined  such  that  all  of  the  physically  relevant  orbits  are  contained 
within  the  interior  of  this  region.  We  define  this  region  as  the  admissible  region ,  and  each  discretized  point  on  the 
admissible  region  as  virtual  debris  particles  (VD’s).  A  set  of  criteria  C  defining  the  admissible  region  has  been 
proposed  by  Tommei,  et  al,  and  later  refined  by  Maruskin,  et  al.  [3]  [6].  This  set  assumes  radar  observations  for 
objects  in  low  and  medium  Earth  orbits  (LEO,  MEO): 


C  = 


(3) 


and 

Ci  =  {(p,p)  :  E  <  0}  C2-  {(p,p)  :  Pmin  ^  P  ^  Pmax) 

C3  =  Up,p)  :  1.03  <  rp)  C4  =  \(p,p) :  ra  <  25},  (4 

where  E  is  the  specific  geocentric  energy  of  the  debris  particle,  pMiN  and  Pmax  are  bounds  for  the  physical  range  of 
the  particle  chosen  “a  priori,”  and  ra  and  rp  are  the  apoapsis  and  periapsis  radii  of  the  orbit  in  units  of  Earth  radii, 
respectively.  In  this  paper,  we  set  (pmin>  Pmax)=(0.3,20)  Earth  radii  as  we  are  interested  in  all  objects  observable  by 
optical  sensors  but  outside  of  the  range  of  radar  sensors,  corresponding  to  an  altitude  of  2000  kilometers  (0.3  Earth 
radii)  to  130,000  kilometers  (20  Earth  radii).  Fig.  1  is  an  example  of  an  admissible  region. 

2.3  Transformation  of  VD’s 

Let  us  consider  the  transformation  of  VD’s  from  topocentric  spherical  coordinates  (i.e.  a,6,a,6,p,p)  into  Poincare 
variables.  Poincare  variables  are  the  non-singular  canonical  counterpart  to  the  equinoctial  orbit  elements  [7].  Their 
main  advantage  is  that  the  variables  can  be  naturally  grouped  into  coordinate-momenta  symplectic  pairs. 
Furthermore,  they  are  defined  and  nonsingular  even  for  circular  and  zero-inclination  orbits.  The  Poincare  elements 
with  respect  to  the  classical  orbit  elements  are  given  as: 


p 

Fig.  1.  An  admissible  region  for  X  =  (a,  6,  a,  6, 0, 0)  =  (2.064,  -0.2378,  0.5072,  0.0654,  0.1,  4.8).  The  different 
shadings  represent  the  different  regions  which  satisfy  each  criterion  in  set  C;  thus,  the  admissible  region  is  where  all 
types  of  shading  overlap,  or  the  region  outlined  by  the  black  line. 


I  =  £l  +  aj  +  M  £  =  yfjld 

g  =  ^2£  (l  -  cos(cj  +  £2)  C5  =  -gtan(U  +  U)  ^ 

l)  =  ^2£  (1  -  cos/)  cos  Q  $  =  -I)  tan  Q, 

A  detailed  discussion  is  given  in  Fujimoto  and  Scheeres  [8]. 

2.4  The  Geometry  of  Admissible  Region  Maps 

Since  the  transformation  from  topocentric  range  /  range-rate  to  Poincare  element  space  is  one-to-one  and  invertible, 
the  map  of  the  admissible  region  is  a  2-dimensional  bounded  submanifold,  or  a  disk,  in  6-dimensional  Poincare 
space  [3]  [8].  Suppose  we  have  maps  from  two  observations  which  have  been  dynamically  evolved  or  regressed  to  a 
common  epoch  r:  FTx{  and  F\2-  From  the  theory  of  general  position,  the  dimension  of  intersection  d  between  a  pair 
of  disks  of  dimension  k  and  /  in  ^-dimensional  space  is  given  as: 


d  =  (k  +  1)  -  n. 


(6) 


where  ifd<0,  the  two  disks  do  not  intersect  randomly  [2].  For  our  problem,  d  =  (2  +  2)  -  6  =  -2,  so  if  two 
admissible  region  maps  intersect  at  all,  it  is  extremely  unlikely  that  they  are  two  separate  objects.  In  addition, 
equation  (6)  indicates  that  if  we  are  to  embed  these  disks  in  5 -dimensional  Poincare  space  such  as  the  (£,  C5,  g,  §>,  b) 
space,  they  will  still  not  intersect  randomly.  Similarly,  for  4-dimensional  Poincare  space  they  intersect  at  a  point, 
for  3 -dimensions  a  line,  and  for  2-dimensions  an  area. 

3.  THE  CORRELATION  OF  OBSERVATIONS  AND  PROBABILITY  DISTRIBUTIONS 

In  this  section,  we  outline  how  to  use  probability  density  functions  of  multiple  data  sets  that  characterize  the  debris 
population  and  combine  them  in  some  standard  comparison  space.  Usually,  this  comparison  space  is  the  Poincare 
element  space  described  in  Section  2.3  or  some  space  derived  from  it.  A  direct  application  of  this  process  is 
determining  whether  a  number  of  observations  are  of  the  same  object,  and  if  they  are,  what  the  approximate  orbital 
characteristics  are. 

Ultimately,  we  would  like  to  know  the  probability  of  an  optically  observed  object,  characterized  by  attributable 
vector  X ,  being  in  the  vicinity  of  some  coordinate  X  in  a  standard  comparison  space  such  as  the  Poincare  orbit 


element  space.  This  probability  can  be  calculated  for  various  values  of  X  across  the  comparison  space,  so  it  is 
beneficial  to  think  of  it  as  an  element  of  a  probability  distribution  function  (pdf).  We  cannot  rationally  determine  to 
any  level  of  useful  accuracy,  however,  neither  this  probability  nor  its  distribution  based  on  one  optical  observation, 
since  X  only  contains  4  variables  (a,8,a,8)  regarding  the  observed  object’s  position  and  velocity,  whereas  6  are 
required  to  fully  describe  the  object's  orbit.  Therefore,  it  is  necessary  to  combine  multiple  observations  of  the  object. 
Here  we  face  another  problem.  In  a  real-world  setting,  there  is  no  guarantee  that  two  or  any  number  of  arbitrary 
observations  are  of  the  same  object;  that  is,  that  they  are  correlated.  We  would  then  like  to  know  if  some  incoming 
data,  such  as  a  new  observation,  is  related  with  the  aforementioned  pdf,  and  if  so,  how  it  affects  the  pdf. 

In  this  paper,  we  consider  the  following  data  sets  regarding  Earth-orbiting  objects  and  their  observations: 

Set  S\\  Past  observation  data  and  debris  distribution  models.  The  United  States  Air  Force  Space  Command 
(AFSPC)  compiles  and  publicly  distributes  Two  Fine  Element  sets  (TFEs)  for  all  known  objects  that  are  in  Earth 
orbit  [1]  [5].  Currently,  the  catalog  consists  of  approximately  14,000  objects.  Fig.  2  shows  a  scatter  plot  of  the 
objects  in  orbit  over  a-e  and  a-i  space. 


Fig.  2.  Distribution  of  known  objects  in  Earth  orbit  over  a-e  (left)  and  a-i  (right)  space. 


Debris  models  incorporate  computer  simulations  to  account  for  objects  too  small  to  be  observed  (generally  smaller 
than  1  meter  for  optical  observations).  MASTER-2005  by  ESA  is  an  example  of  such  a  model.  The  data  in  this 
category  are  discrete,  and  are  independent  of  X. 


Set  S2(X):  The  distribution  of  virtual  debris  (VD)  particles  over  the  standard  comparison  space.  Although  VD's  are, 
by  definition,  uniformly  distributed  in  the  admissible  region  (i.e.  range  /  range-rate  space),  for  common  choices  for 
the  comparison  space  such  as  orbit  element  and  Poincare  space,  the  distribution  of  VD's  in  such  spaces  is  non- 
uniform  due  to  the  non-linearity  of  the  mapping.  As  a  consequence,  certain  values  of  X  become  more  likely  than 
others.  The  data  in  this  category  are  continuous,  and  are  a  function  of  X.  Computationally,  however,  a  large  and 
discrete  sample  set  (§2)  is  used  instead. 

We  discretize  the  standard  comparison  space  into  M(  =  n6_,  M; )  6-dimensional  hypercubes  (or  “bins”),  which  we 

J—L  J 

index  with  vector  i.  By  doing  so,  we  rationally  group  TLE  objects,  modeled  debris,  and  VD's  with  similar  orbital 
characteristics.  In  the  comparison  space,  objects  in  a  particular  bin  are  spatially  indistinguishable;  i.e.  we  treat  their 
coordinates  as  being  the  same  as  those  that  define  the  position  of  the  bin.  The  discretization  makes  up  for 
deficiencies  in  5)  data  and  undersampling  of  S2,  as  well  as  speed  up  computational  turnaround.  Note  that  Mis  an 
important  parameter,  as  if  it  is  too  small,  we  lose  too  much  of  the  spacial  resolution  of  the  data  sets.  If  Mis  too  large, 
then  the  aforementioned  data  deficiencies  and  undersampling  can  negatively  influence  our  results. 


With  this  discretization,  it  is  natural  to  consider  the  data  in  sets  5)  and  S2  as  discrete  pdfs  spanning  the  comparison 
space  rather  than  a  set  of  countable  elements.  We  refer  to  these  pdfs  as  Si(i)  and  s2(i,  X ),  respectively.  Practical 
definitions  for  and  s2  are  given  by  first  defining  the  following  sets: 


A[  =  {a  :  a  e  S  \  and  is  mapped  to  bin  i} 
#i(£)  =  {b  :  b  e  S 2 (£)  and  is  mapped  to  bin  i). 


(7) 


Then, 


s  i(i)  = 


n(AQ 
n(S  i) 


J2(i,3£) 


n(*i(*)) 

*CS2(*))’ 


(8) 


where  n(A i)  is  the  number  of  elements  in  set  and  so  on. 

Suppose  we  have  some  prior  discrete  pdf  g0(i,r)  that  describes  the  probability  that  a  particular  object  of  interest  O  is 
consistent  with  bin  i  at  epoch  r.  That  is,  if  P  is  a  probability  measure  and  E0( i,  r)  is  an  event  where  O  is  consistent 
with  i  at  time  r,  then: 


P[Eo(lr)\=go(lT).  (9) 

g  can  originate  from  the  TLE  catalog  or  a  debris  distribution  model  (Si),  one  observation  or  a  set  of  observations  that 
we  believe  are  correlated  a  priori  with  O  (S2),  or  any  combination  of  Si  and  S2.  Here,  we  treat  i  as  a  random  variable 
that  spans  the  bin  index  space.  Also,  all  information  has  been  propagated  to  r. 

Now,  let  us  consider  a  new  series  of  information  {r}  that  has  come  in  from  Si  or  S2,  such  as  uncorrelated 
observations.  We  would  like  to  calculate  a  posterior  pdf  h0( i,  r)  on  whether  O  is  consistent  with  bin  i  and  is  related 
to  {r} .  It  is  obvious  that  h0(i,  r)  =  0  if  the  new  information  does  not  regard  O.  In  order  to  filter  out  such  trivial  cases, 
we  let  event  xf  be  one  where  the  series  {r}  is  related  with  O  and  add  this  event  as  a  condition  to  h.  Using  Bayes’ 
theorem: 


P[E0(i,T)\xf]  =  h0(i,T)  = 


/w(i,T)go(i,T) 
Zj  /(r)Cj>  T)go(j,  t)} 


(10) 


where  the  sum  in  the  denominator  is  over  all  bins,  and /{r}(i,  t)  is  a  pdf  that  describes  the  probability  that  {r}  is 
consistent  with  bin  i: 


P[E{r}(i,T)]=f{r}(i,T). 


(ii) 


E{Ty( i,  r)  is  an  event  where  {r}  is  consistent  with  i  assuming  all  information  has  been  dynamically  evolved  to  r. 
Again,  i  is  treated  as  a  random  variable.  If  the  information  in  {r}  are  observations, /is  the  admissible  region  of  that 
observation  mapped  to  the  comparison  space  and  to  epoch  r.  f=  s2( i,  £)•  Similarly,  if  the  information  is  the  TLE  set, 
then /=  si(i).  Note  that  from  (6),  if/>  0  and  g  >  0  at  some  bin  i  regardless  of  discretization  size,  then  {r}  and  O  are 
most  likely  related.  Therefore,  given  O  is  consistent  with  i,  whether  {r}  and  O  are  related  depends  only  on  whether 
{r}  is  consistent  with  i: 


P[x?\E0(i,T)\=Mi,T).  (12) 

Furthermore,  the  converse  of  the  above  argument  ensures  that  as  long  as  {r}  and  O  are  related,  then  2j  /  ‘  8  >  0,  so 
(10)  is  well-defined. 

In  a  graphical  sense,  pdf  h  is  a  "cut  out”  of  the  region  where / and  g  intersect;  h  >  0  for  any  bins  where  both />  0 
and  g  >  0,  and  the  probability  expressed  by  h  is  one  that  is  evaluated  over  this  overlap  region.  Based  on  (6),  we  can 
look  at  whether  h  >  0  for  some  bin  i  to  deduce  with  confidence  whether  or  not  the  new  information  is  related  to  the 
object  of  interest. 


4.  NUMERICAL  SIMULATIONS 

In  this  section,  we  discuss  results  from  an  implementation  of  our  method  in  MATLAB.  Note  that  our  goal  is  to 
show  that  our  method  is  able  to  correlate  objects  and  give  an  initial  orbit  estimate  by  giving  it  only  angle  and  angle- 
rate  information  from  the  observations.  Determining  its  robustness  in  real-world  situations  is  future  work. 


We  extracted  8  objects  from  the  two-line  element  (TLE)  catalog  to  obtain  a  sample  set,  and  refer  to  them  as  follows 

[1]: 


•  3  objects  in  GEO  (GEOl-3), 

•  1  object  in  a  Molniya  orbit  (MOL1), 

•  2  object  in  an  eccentric  MEO  orbit  (EM1,EM2), 

•  1  object  in  a  circular  MEO  orbit  (CM1),  and 

•  1  GPS  satellite  (GPS1). 

The  orbital  parameters  of  each  object  is  given  in  Appendix  A. 

We  simulated  996  zero-error  observations  of  right  ascension,  declination,  and  their  time  derivatives  made  from  4 
observatories  for  all  8  objects  over  the  course  of  24  hours.  Thus,  £  =  (a,  8 ,  a,  8,  to,  0, 0),  where  t0  is  the  observation 
epoch  and  (0,  cp)  is  the  inertial  angular  position  of  the  observation  point.  We  can  assume  zero-error  in  the 
observations  because  the  uncertainty  in  the  angular  information  is  generally  much  less  than  the  uncertainty  in  the 
range  and  range-rate;  refer  to  Fujimoto  and  Scheeres  for  a  discussion  on  the  effects  of  observation  error  to  the 
outcome  of  our  method  [8].  The  4  observations  points  were  located  at: 

•  Socorro,  New  Mexico  (33.8172°N  106.6599°W) 

•  Maui,  Hawaii  (20.7088°N  156.2578°W) 

•  Diego  Garcia,  British  Indian  Ocean  Territory  (7.41 173°S  72.45222°E) 

•  Moron  Air  Base,  Spain  (37.170°N  5.609°W) 

Approximately  every  15  minutes,  the  code  generates  attributable  vectors  (i.e.  simulated  observations)  for  all  objects 
that  are  above  the  local  horizon  at  any  given  observation  point. 

We  assumed  no  a  priori  information  regarding  the  observed  objects,  and  thus  used  a  uniform  initial  pdf.  The 
discretization  of  the  Poincare  space  was  such  that: 
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and  (100,77,123,123,123,123)  bins  in  each  coordinate  direction  for  a  total  of  1.7624><1012  bins.  All  dynamics  were 
two -body. 

Fig.  3  is  a  graphical  representation  of  the  process  explained  in  Section  3  for  2  observations.  The  red  and  green 
regions  each  represent  pdfs  based  on  observations  that  have  been  dynamically  evolved  to  a  common  epoch  (i.e. 
pdfs  /  and  g).  The  propagation  has  “shredded”  the  red  pdf  in  the  £-1  plane  [3].  The  blue  region  is  the  combined 
distribution  (i.e.  pdf  h).  The  yellow  asterisk  is  the  true  state  of  the  observed  object.  The  distributions  have  been 
projected  onto  2-dimensional  subspaces  using  their  coordinate-conjugate  momentum  pairs;  note,  however,  that  the 
correlation  was  conducted  in  the  full  6-dimensional  Poicare  space.  When  correlating  two  observations  of  the  same 
object  (top),  we  see  that  h  >  0  for  a  very  small  region  of  the  state  space;  for  this  particular  example,  h  >  0  for  1 1 
bins.  Furthermore,  the  true  state  is  included  in  the  region  in  state  space  where  h  >  0.  Therefore,  the  state  estimate  is 
good.  On  the  other  hand,  when  two  observations  are  of  different  objects  (bottom),  h  =  0  for  the  entire  state  space, 
which  allows  us  to  conclude  that  the  two  observations  are  unrelated. 


^  ^ 


Fig.  3.  Projections  of  probability  distributions  when  correlating  observations  of  the  same  object  (top;  EMI -EMI) 
and  a  different  object  (bottom;  EM1-EM2).  Length  units  in  Earth  radii,  time  units  in  hours,  mass  units  in  kilograms, 

angle  units  in  radians. 


Our  correlation  process  performed  well  for  all  996  observations:  all  observations  were  correctly  correlated  to  the  8 
objects  and  their  states  were  correctly  estimated  down  to  a  region  of  at  most  4  bins.  Thus,  there  were  no  “false 
positive”  results.  We  developed  a  linearization  technique  in  which  the  pdfs  are  mapped  linearly,  reducing  the 
computational  time  while  attaining  good  accuracy  [8].  On  a  dual-core  Xeon  server  with  32-bit  numerics,  each 
correlation  run  took  approximately  10  minutes.  If  we  wish  to  further  reduce  the  region  over  h  >  0  as  well  as  reduce 
computation  time,  we  can  assume  a  priori  that  all  observed  objects  were  included  in  either  the  TLE  catalog  or  some 
debris  distribution  model  instead  of  the  uniform  distribution  assumption  we  made  for  Fig.  3.  Then,  the  admissible 
region  maps  are  “pre-conditioned”  to  exclude  unrealistic  objects.  Correlation  times  were  reduced  to  1  to  2  minutes. 
Table  1  lists  the  number  of  overlap  bins  for  each  correlated  object. 
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MOL1 

1 
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1 

YES 

EMI 

2 
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1 
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EM2 

3 
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1 
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CM1 

2 
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1 
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GPS1 

3 

YES 

1 

YES 

Table  1.  The  estimation  accuracy  of  each  observed  object  when  using  a  uniform  initial  pdf  (left)  and  the  TLE  as  the 

initial  pdf  (right). 

5.  LIMITING  CASES 

In  this  section,  we  investigate  special  observations  cases  where  we  expect  the  method  to  have  difficulty  in 
calculating  an  accurate  initial  orbit  estimate;  namely,  when  two  objects  lie  within  a  bin  size,  when  an  object  is 
observed  simultaneously  at  two  observation  points,  and  when  objects  in  a  satellite  constellation  are  observed  over 
long  periods  of  time. 


5.1  Limitation  of  Bin  Resolution 

Recall  from  Section  3  that  all  objects  within  a  bin  are  regarded  as  having  the  same  state  parameters.  Therefore, 
there  may  be  cases  where  two  satellites  that  are  close  in  state  space  are  falsely  correlated.  We  can  always  overcome 
this  limitation  by  increasing  the  fidelity  of  the  discretization;  one  efficient  way  is  to  implement  a  recursive  algorithm 
where  the  bin  size  is  decreased  in  the  vicinity  of  overlapping  region.  Moreover,  we  can  be  smart  about  how  we 
initially  discretize  the  space.  For  instance,  we  expect  that  a  finer  grid  is  necessary  in  the  I  direction  (i.e.  mean 
anomaly)  compared  to  the  £  direction  (i.e.  semi-major  axis). 

5.2  Simultaneous  Observations  of  an  Object 

The  proposed  method  may  run  into  difficulties  computing  a  precise  state  estimate  when  the  two  pdfs  are  near 
parallel.  Although  such  pdfs  would  still  most  likely  intersect  at  a  single  point,  they  may  appear  to  occupy  the  same 
bins  in  the  discretized  state  space  (i.e.  overlap)  over  a  large  region.  One  case  where  pdfs  become  near  parallel  is 
when  an  object  is  observed  simultaneously  at  2  different  observations  points,  as  seen  in  Fig.  4.  Here,  pdf  h  spans 
over  688  bins  using  the  nominal  discretization  in  Section  4.  This  result,  however,  does  not  imply  that  the  two  pdfs 
intersect  over  a  large  planar  region;  from  (6),  2-dimensional  intersections  of  pdfs  are  extremely  unlikely.  Indeed,  as 
we  refine  the  discretization  by  20%,  then  50%,  the  overlap  region  begins  to  converge  upon  the  coordinate  of  the  true 
object  state.  Note  that  compared  to  the  method  proposed  by  Maruskin,  et  al.  which  evaluates  intersections  of 
manifolds  within  their  2-dimensional  projections,  the  new  method  converges  faster  as  the  intersections  are  evaluated 
in  the  full  6-dimensional  space  [3]. 


Fig.  4.  Combined  pdf  h  of  two  near-parallel  pdfs  (CM1  observed  simultaneously  at  (@,<p)  =  (1,1)  and  (1.5,2)  at  time 
t  =  26.36  hours)  as  the  discretization  is  refined  from  M=  1.7624x10 12  (top)  to  1.2M  (middle)  and  1.5M  (bottom). 


5.3  Observations  of  Satellite  Constellations 

Suppose  we  make  two  observations  separated  by,  say,  4  hours,  which  are  consistent  to  the  same  orbit  plane.  We  can 
either  conclude  that  we  observed  1  object  with  an  orbital  period  that  is  any  divisor  of  4,  or  that  we  observed  2 
different  objects  in  a  satellite  constellation  with  an  orbit  period  that  is  any  divisor  or  multiple  of  4.  For  this 
particular  example,  most  of  the  solutions  in  the  former  set  will  have  a  semi-major  axes  that  are  too  small  to  be 
included  in  the  admissible  region.  If  the  temporal  separation  of  the  observations  were  larger,  however,  then  single¬ 
object  solutions  may  become  viable.  The  long-term  propagation  “dilutes”  information  regarding  the  object’s  angular 
position  that  we  can  extract  from  the  observation  separation  time. 

Mathematically,  let  two  object  be  on  an  orbit  with  period  r  separated  by  mean  anomaly  AM.  In  two-body  motion,  r 
is  related  to  semi-major  axis  as  r  =  2 n^a3/n.  Suppose  we  observe  one  of  the  objects  at  time  0,  and  then  the  other  at 
time  Nt  +  ( AM/n ),  where  N  £  N  and  n  is  the  mean  motion;  i.e.  we  observe  the  second  object  after  N revolutions. 
Now,  if  we  were  to  wrongly  assume  that  we  observed  the  same  object  twice,  then  the  true  anomaly  a  +  Aa  of  this 
ficticious  object  is: 


Therefore,  N  — »  oo  =>  Aa/a  —>  0.  If  the  first  observation  generated  a  non-empty  admissible  region,  then  it  is  likely 
that  the  proposed  method  will  mistakenly  relate  the  two  observations  given  they  are  temporally  well-separated. 


Fig.  5  is  a  graphical  representation  of  the  above  scenario.  We  simulated  observing  2  separate  objects  in  a 
constellation  from  one  observatory  with  different  observation  separation  times.  When  the  separation  is  14.85  hours, 
the  combined  pdf  is  null  at  all  bins,  meaning  the  method  successfully  recognizes  the  observations  as  that  of  different 
objects.  When  the  separation  is  increased  to  60.61  hours,  however,  we  see  that  the  2  pdfs  overlap,  and  thus  we 
obtain  a  “false  positive”  result. 


Fig.  5.  Correlation  of  observations  of  2  different  objects  separated  by  a  mean  anomaly  of  2/371  on  the  same  orbital 
plane  as  EM2.  Observations  are  separated  by  14.85  hours  (top)  and  60.61  hours  (bottom).  All  observations  are 

made  at  (0,  (p)  =  (2,4). 


6.  CONCLUSIONS 


In  this  paper,  we  discussed  methods  of  correlating  multiple  optical  observations  as  well  as  providing  initial  state 
estimates  using  pdfs  in  the  Poincare  orbit  element  space.  We  outlined  the  method  that  incorporates  Bayes’  rule.  An 
implementation  of  the  method  in  MATLAB  successfully  correlated  996  optical  observations  and  provided  good 
initial  orbit  estimates.  Finally,  some  limiting  cases  were  examined. 

Future  work  will  be  incorporate  more  accurate  dynamics,  and  to  apply  our  method  to  different  observation  scenarios 
such  as  space-based  observations  and  observations  of  objects  in  heliocentric  orbit. 

7.  APPENDIX  A:  ORBITAL  ELEMENTS  FOR  OBJECTS  IN  THE  NUMERICAL  EXAMPLE 

The  classical  and  Poincare  orbital  elements  of  all  8  objects  in  our  simulation  sample  set  are  listed  below  as: 

(a  [rE],  e,  i  [rad],  Q  [rad],  [rad],  M  [rad]),  (fi  [rE2/hour],  I  [rad],  ©  [rE/hour1/2],  g  [rE/hour1/2],  §  [rE/hour1/2],  \)  [rE/hour1/2]), 
where  rE  is  Earth  radius. 

•  GEO 

GEOl  (6.6102,0.0003,0.0002,3.1274,2.3294,5.7226)  (11.4721,4.8962,0.0006,0.0006,-0.0000,-0.0006) 

GE02  (6.6109, 0.0003, 0.0009, 5.9501, 2.9681, 5.8729),  (11.4727, 2.2247, -0.0005, -0.0009, 0.0010, 0.0029) 

GE03  (6.6109,0.0001,0.0001,3.0902,4.2464,6.1 176),  (11.4727,  0.8879,-0.0003,0.0001,-0.0000,-0.0002) 

•  Molniya 

MOL1  (4.1971,0.7154,1.1204,0.8126,5.1137,0.1671),  (9.1414,6.0934,0.8200,2.1991,-1.9501,1.8469) 

•  Eccentric  MEO 

EMI  (3.6573,0.7123,0.3106,1.0976,1.6080,5.9942),  (8.5333,2.4166,-0.9526,-2.0447,-0.6739,0.3451) 

EM2  (4.1472,0.5529,1.2347,5.5811,2.4137,4.8996),  (9.0868,0.3281,-1.7237,-0.2444,2.0573,2.4323) 

•  Circular  MEO 

CM1  (3.9994,0.0006,1.1284,4.9148,4.2128,2.9461),  (8.9234,5.7905,-0.0005,-0.0017,3.1297,0.6421) 

•  GPS  satellite 

GPS1  (4.1645,0.0048,0.9599,2.7242,3.6934,2.5851),  (9.1057,2.7195,-0.0019,0.0143,-1.1296,-2.5474) 
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